Ordering And Partial Ordering In Holmium Titanate And Related Systems 
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Abstract 



Wc take another look at two compounds which have been 
discussed as possible realizations of "spin ice", namely 
holmium titanate and dysprosium titanate. As we have 
earlier observed, holmium titanate does not display ice- 
like behaviour at low temperatures because the long ranged 
dipolar interactions between spins are strong compared to 
the nearest neighbour interactions. We show, analytically, 
that the true ground state of this system must be fully or- 
dered, but simulations only reach partially ordered states 
because there are infinite energy barriers separating these 
from the true ground state. We also show that the true 
ground state of our model of dysprosium titanate is also 
fully ordered, and offer some explanations as to why sim- 
ulations and experiments show ice-like behaviour. We dis- 
cuss the effect on these systems of an applied magnetic 
field. Finally, we discuss several other models which show 
similar partial or full ordering in their ground states, in- 
cluding the well known Ising model on the fee lattice. 

1 Introduction 

With the wide interest in the physics of disorder and 
frustrated magnetism, pyrochlore magnets have attracted 
great attention in recent years [jj, and it is especially 
interesting to consider pyrochlores well approximated by 
the fsing model. We recently discovered [|| that dys- 
prosium titanate, an fsing pyrochlore, exhibits a ground 
state entropy very much like that of ice. Anderson had 
long ago predicted that this should happen for a nearest- 
neighbour fsing pyrochlore. However, the story is not quite 
so simple here: the dominant interaction is really a long- 
ranged dipole-dipole interaction. Moreover, the similar 
Ising pyrochlore holmium titanate has very different low- 
temperature behaviour. We had explored the reasons for 
this in our earlier papers 0, pf; here we take those argu- 
ments further. 

We observed earlier M that simulations of a model of 
the Ising pyrochlore holmium titanate suggest that it has 
a partially ordered ground state, which is degenerate but 
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has no entropy per particle (the degeneracy being of the or- 
der 2 L where L is the system length, rather than 2 L ), and 
there appears to be a first order phase transition from the 
high temperature paramagnetic phase to the above par- 
tially ordered phase. In this article we clarify the nature 
of this partial ordering, and explain it without recourse 
to simulations. Actually, we show that the ordering of 
the true ground state here is in fact complete, but there 
are numerous low-lying partially ordered metastable states 
which are separated from the true ground state by infinite 
energy barriers. It is easy for the simulation (and, pre- 
sumably, the real compound) to get stuck in one of these 
states on cooling, and impossible to reach the true ground 
state in finite time thenceforth. Moreover, this is also the 
true ground state of our model of dysprosium titanate; but 
both the model and the experiments [pi suggest a finite 
ground state entropy characteristic of nearest-neighbour 
"spin ice". This, we suggest, is because this system has 
stronger nearest-neighbour interactions and settles into an 
ice-like state at a higher temperature (over 1 K), and can- 
not then be easily dislodged from this into the true ground 
state. Our exact results, and also our simulations, support 
and substantiate our original suggestions of a transition to 
partial ordering, in contrast to recent suggestions to the 
contrary |4(], namely that that our model for holmium ti- 
tanate should have an ice-like ground state. 

We also substantiate the major results and the underly- 
ing model from our earlier work [Bj , namely that there is a 
low temperature entropy observed in dysprosium titanate, 
it is decreased in the presence of a magnetic field, and the 
interactions in the system are dipole-dipole magnetic in- 
teractions and an isotropic superexchange. The fact that 
our simulations with a magnetic field reproduce experi- 
mental results quite well, qualitatively and quantitatively, 
confirms our calculation of the dipole moment (which is 
the only thing that couples with the field) and our esti- 
mate of the superexchange. Moreover, the ground state in 
the presence of a strong field is not the same as the zero- 
field ground state. While the experiments were done on 
powdered samples, simulations suggest that the behaviour 
of dysprosium titanate in a magnetic field is very direction 
dependent. A strongly direction dependent ordering was 
initially suggested on the basis of the specific heat mea- 
surement done on a powdered sample H. We compare 
the experimental data on powder samples with simulations 



averaged over large numbers of directions, compare sharp 
features in both, and make important predictions for pos- 
sible future experiments on single crystals. 

Next, we use the insight from the ground state analysis 
to progressively reduce the pyrochlore Ising model to a se- 
quence of simpler models which display similar behaviour: 
namely, a six-vertex model on the "diamond lattice" with 
non-local interactions, which reproduces all the essential 
behaviour of holmium titanate, and has a fully ordered 
ground state and several metastable partially ordered low- 
lying states (section || a six-state magnetic model on the 
fee lattice, which actually has the sort of partially or- 
dered ground state that the simulations had suggested 
for holmium titanate, and also reproduces the important 
behaviour of holmium titanate (section H); and the well- 
known Ising model on the fee lattice, which has been stud- 
ied before [pf and is known to have exactly the same sort 
of partial ground state ordering that concerns us here (sec- 
tion |7|) . Along the way, we also introduce a square-lattice 
vertex model, by analogy with the above diamond-lattice 
model, which may be worthwhile to study in its own right. 

Actually, the simplest example of partial ordering in an 
Ising system is perhaps the triangular lattice antiferromag- 
net, with interaction Jj along bonds in one direction (say 
parallel to the x axis), J2 < J\ along bonds in other direc- 
tions (Figure U\) , and J\ , J2 positive (antifcrromagnetic) . 
Then we have perfect antiferromagnetic order along a line 
of sites in the x direction, but each atom on an adjacent 
chain is frustrated so the adjacent chain (which also has 
perfectly antiferromagnetic order) has two possible con- 
figurations with respect to the first, and the system as 
a whole has a degeneracy 2 L where L is the number of 
chains. This system is exactly solvable M, and has no 
finite-temperature phase transition. This situation is quite 
relevant to what happens in our system. The analogous 
model in three dimensions, the fee Ising model, is discussed 
at the end of the paper and in the cited references. 
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Figure 1: Anisotropic triangular lattice, with antiferro- 
magnetic Ising interaction J along solid lines and J' < J 
along dashed lines. In ground state, each chain along solid 
lines is ordered, but there is no order along dashed lines. 



2 Ordering in our model of 
Ho 2 Ti 2 7 

We briefly recapitulate our model of holmium titanate. 
The underlying lattice is the pyrochlore lattice, a lattice 
of corner-sharing tetrahedra of two possible orientations, 
which is well visualised as an fee lattice of tetrahedra (Fig- 
ure 0). It may be generated by taking a single tetrahe- 
dron of one orientation and translating it by the primi- 
tive basis vectors of the fee lattice (Figure 0); the tetrahe- 
dra of the other orientation emerge automatically by this 
procedure — see Figure g. Thus we use the lattice vectors 
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Figure 2: The conventional unit cell of the pyrochlore, a 
lattice of corner-sharing tetrahedra with fee symmetry. 

We have Ising spins (the / electron states of the holmium 
atoms) located at these points. The local Ising axis is the 
line joining the centres of the two adjoining tetrahedra: 
thus, each spin points directly out of one tetrahedron, and 
into the next one. The spins carry magnetic moments cor- 
responding to J = 8, g s (the Lande factor) = 1.25. Based 
on this, the expected nearest-neighbour dipole-dipole in- 
teraction has an energy of ±2.35 K, the sign depending on 



the alignment of the spins: one pointing out, one in is pre- 
ferred. However, the experimental compound (unlike its 
cousin, dysprosium titanate) has properties which we can 
only explain by postulating a significant nearest-neighbour 
superexchange, which we estimate from high temperature 
expansions for the susceptibility to be around 1.9 K with 
opposite sign to the dipolar interaction. So, with the con- 
vention that an Ising spin S = 1 if it points out of an "up" 
tetrahedron and S = — 1 otherwise, we write the Hamilto- 
nian as follows: 
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0.45 Kelvin, for nearest neighbour spins (4) 
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where n, is a unit vector pointing along the Ising axis 
at site i in the outward direction from an "up" tetrahe- 
dron. Thus, this system has a drastically reduced nearest- 
neighbour interaction energy, which means the importance 
of the further-neighbour interactions increases. We saw 
that with only a long ranged dipole-dipole interaction be- 
tween these Ising spins, the behaviour seems to change 
little from the nearest-neighbour Ising model which has a 
finite ground state entropy; we speculate that the substan- 
tially different behaviour of holmium titanate is due to the 
significantly greater importance of further-neighbour in- 
teractions. As noted earlier, simulations of such a model 
do predict a freezing of the system at around 0.7 K, in 
agreement with experiment. 




Figure 3: The pyrochlore lattice can be generated by 
translating a tetrahedron along vectors ai, a2, a3. Tetra- 
hedra of the opposite orientation are formed from the cor- 
ners of these (thin black lines). 

We now find the ground state (GS) for this system. For 
clarity we refer to the two different orientations of tetra- 
hedra that occur as "up" and "down" : each "up" tetrahe- 
dron shares corners with only "down" tetrahedra, and vice 
versa. Consider a single tetrahedron: it has six possible 



ground states, each of which has two spins pointing out and 
two in. We label these states A, A', B, B', C, C where A' 
is A with the directions of all spins reversed, and similarly 
B' and C (Figure 0) . Consider a cluster of sites consist- 
ing of a single "up" tetrahedron and its translation by the 
three primitive lattice vectors of the fee lattice (as in Fig- 
ure ||) , and long-ranged dipole-dipole interactions spanning 
the entire cluster plus a nearest-neighbour superexchange, 
as described above. We find the GS of this "tetrahedron 
of tetrahedra" by the unobjectionable method of enumer- 
ating each of the 2 16 allowed states on a computer, and 
picking out the lowest-energy ones. 
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Figure 4: Three possible ground state configurations A,B, 
C of a tetrahedron. We will denote by A', B', C" these 
respective configurations with all spins reversed. '+' in- 
dicates a spin pointing out of the tetrahedron, '—'a spin 
pointing into the tetrahedron. 
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Figure 5: Two of the twelve allowed ground states for 
a cluster of four tetrahedra. Each of the twelve states is 
related to the others by rotations or reflections. This con- 
sideration by itself leads to partial ordering in the ground 
state of the full system, as described in the text. 

We find that the GS of this cluster is twelve-fold de- 
generate. Two are shown in Figure 0. All others are ro- 
tations/reflections of these (which are reflections of each 
other). In each GS, as one might expect, only states sat- 
isfying the ice rule occur. Furthermore, only two kinds of 
ice-like configurations occur: A and A' , or B and B' , or C 
and C" . They occur twice each. In the case shown in the 
figure, once the configuration of tetrahedron at the origin 
is fixed as C, the tetrahedron at a3 must have the op- 
posite configuration (C); and the tetrahedron at a 2 may 
have either configuration (C or C"), but tetrahedron a! 
must have the opposite configuration to a 2 . With the pre- 
viously quoted values for the dipole and superexchange 
interactions, the configurations in Figure have an energy 



—7.5 K, and the next lower energy configurations have an 
energy —6.9 K. 

If we note that each tetrahedron in its ground state has 
a dipole moment, which is perpendicular to two possible 
lattice translation vectors, the rule is this: in every GS 
only configuration of two opposite kinds occur (say C and 
C"), which therefore have antiparallcl magnetic moments, 
and two tetrahedra separated by a vector perpendicular to 
these moments must have opposite configurations. These 
are the twelve ground states allowed for this cluster, but we 
are making no theoretical argument for this; this is what 
we learn from brute-force enumeration of states. In any of 
these states, moreover, the included "down" tetrahedron 
also has an ice-ruled configuration. 

Now consider the entire fee lattice of tetrahedra (that 
is, the pyrochlore lattice). Once we fix the configuration 
of a single "up" tetrahedron (say, as C) each of its neigh- 
bouring tetrahedra, with which it forms part of a cluster 
of the above sort, must have the same or opposite con- 
figuration (say C or C); and by extending further from 
each of these, this must be true for the entire lattice. We 
can also see that if we travel in either of the two lattice 
directions perpendicular to the magnetic moment of these 
tetrahedron configurations, we must have a perfectly alter- 
nating sequence (say C, C", C, C", . . . ); but when travel- 
ling in a third direction we have no ordering rule. So we 
get a ground state ordering in two directions but not in 
the third, and a degeneracy exponential in L (the system 
length) rather than L 3 (the system volume), precisely as 
the simulations suggested. 

This is not the whole story, though. The system can be 
equally well described in terms of "down" tetrahedra, so 
the same sort of ordering should be evident if we describe 
a ground state configuration using "down" tetrahedra. So 
only two configurations for these tetrahedra should be al- 
lowed, each of which is the other with all spins reversed. 
But we note immediately that the two "down" tetrahedra 
displayed as grey lines in the two clusters in Figure || do 
not have opposite configurations. It follows that the two 
cluster configurations in that figure cannot both occur in 
the ground state: only one can. This immediately implies 
that the ordering sequence in a third direction is not ran- 
dom, but constant (say, C, C, C, . ..). 

So the true ground state for our model of holmium ti- 
tanate is only twelvefold degenerate, and viewed in terms 
of configurations of either upward or downward tetrahedra, 
consists of alternating ordering of opposing configurations 
in two directions but a constant configuration in a third di- 
rection. However, the partially ordered states are also very 
low in energy, and moreover a system stuck in such a state 
can only get out by flipping entire planes of tetrahedron 
configurations, which is impossible in the thermodynamic 
limit. So simulations tend to get stuck in such states and 
in other "domainised" states and the chances of a given 
simulation actually hitting a true ground state are very 
small — unless we use some sort of specialised "cluster" al- 
gorithm which may not imitate the dynamics of the real 



system very well. 

This is exact except for one thing: we have ignored in- 
teractions between further-neighbour tetrahedra, so effec- 
tively confined our interaction range to the 5th neighbour, 
which is the maximum separation of spins in two adja- 
cent tetrahedra. As long as only nearest-neighbour "up" 
tetrahedra interact, so that the whole system really can 
be decomposed into clusters as in Figure g, our picture 
is totally correct. In the presence of long ranged interac- 
tions, luckily the nature of the system (with four different 
local z directions) is such that the effects of the more dis- 
tant tetrahedra cancel heavily and have little effect on the 
energy of a single cluster. Thus, if one take a particular 
"up" tetrahedron in the ground state, its energy turns out 
to be —21 K because of interactions with the immediately 
neighbouring "up" tetrahedra (to which the cluster argu- 
ment applies), but only —0.12 K because of interactions 
with all other tetrahedra in the system. The effect is not 
only small, but tends to stabilise this order. (With a ran- 
dom ice-like configuration, this energy averages to zero, but 
fluctuates considerably from site to site.) Recall further- 
more that the cost of disturbing a single four-tetrahedron 
cluster from its ground state (Figure ||) is at least 0.6 K, 
and in fact much more since each tetrahedron is shared by 
four such clusters. 

So we can be satisfied that the long ranged interactions 
will not affect the above conclusions. In fact, the simula- 
tions too don't show much dependence on the range of the 
interaction, provided it extends to at least the third neigh- 
bour (Figure |6|) , as indeed we had argued in our earlier pa- 
per, where we cut off the interaction at the fifth-neighbour 
distance |2|. The suggestion [Q that ice-like behaviour is 
restored by cutting off after the 10th neighbour seems un- 
tenable to us, and we do not observe it in our simulations 
even on extending the interaction to the 12th neighbour 
(which is halfway across our sample) . Possibly the different 
results obtained in H is due to additional approximations 
involved, such as Ewald sums, and a somewhat smaller 
sample size; our simulations use no approximations in cal- 
culating the energy except the cutoff, which as we have 
seen, is quite justifiable. Longer simulations on bigger sys- 
tems may throw more light on this question, but we now 
turn to some other interesting aspects of the problem. 

This true ground state ordering is more easily visualised 
(though less easily analysed) with the cubic unit cell rather 
than our parallelepiped; this is discussed in section |[ 

3 Dysprosium titanate 

Dysprosium titanate, which we earlier reported as show- 
ing ice-like behaviour experimentally and in simulations, 
appears to have a much weaker superexchange between 
nearest neighbours. With a nearest-neighbour-only model 
of the superexchange, we find we need a superexchange 
of around +1.1 K, that is, the nearest-neighbour interac- 
tion is around —1.25 K compared to the bare dipole-dipole 
value of —2.35 K. With these numbers, we get a reason- 
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Figure 6: The simulated specific heat when the interaction 
is cut off at the third {R < 2r), fifth (R < 2.646r), and 
twelfth (R < 4r) nearest neighbour distances (r = 3.54 A 
roughly). The position of the phase transition and the 
plot of the specific heat near the transition hardly change 
at all on increasing the range of the interaction; but one 
needs longer equilibriation times with increased interaction 
ranges. The simulations show a significant energy drop at 
the transition, suggesting that it is first order. 



able agreement of the simulation with experiment (Fig- 
ure f 7 ]). However, we get even better agreement by uni- 
formly scaling down the long-ranged dipole-dipole interac- 
tion by roughly this factor. (This is what was reported in 
our earlier paper.) This suggests that the superexchange is 
not strictly nearest-neighbour but extends over to further 
neighbours. 

We need to understand why Ho2Ti207-like behaviour is 
not observed in this case. With either of the two mod- 
els above (small nearest-neighbour-only superexchange, 
or uniformly scaled-down dipole-dipole interaction) the 
ground state for the cluster of four tetrahedra remains the 
same as before, the difference in energy from the next- 
lowest state too remains roughly the same, and the above 
arguments should still go through. 

The difference is that, because of the stronger nearest- 
neighbour interactions here, this system undergoes a 
crossover from a paramagnetic phase to an ice-ruled phase 
at a considerably higher temperature (greater than 1 K); 
and by the time it cools down to the temperature (< 0.7 K) 
where we expect a transition of the sort described here, 
it is already stuck in a disordered ice-like state and (be- 
cause of the stronger nearest-neighbour interactions) can- 
not easily break out of this state to access other states. 
It appears that the temperature of the crossover to the 
ice-like phase is dictated by the nearest-neighbour inter- 
actions. In dysprosium titanate these have an energy of 
around 1.3 K, and hence the ice rule is already in place 
by the time we go down to 0.7 K and the spins are al- 
most frozen, thus the ordering transition no longer has a 



Figure 7: A comparison with experimental data for dys- 
prosium titanate of a simulation using nearest-neighbour 
superexchange and long-ranged dipole-dipole interactions, 
and a simulation using uniformly scaled down dipole-dipole 
interactions (as was done in ref. 0]). 



chance to occur. In holmium titanate, on the other hand, 
the nearest-neighbour interaction is around 0.4 K, thus at 
0.7 K the system is in no sense frozen, plenty of spin flips 
take place and the ordering transition occurs (Figure 0). 

So while our arguments show that the true ground state 
here is ordered and only twelve-fold degenerate, the system 
tends to get stuck in fairly generic ice-like states. We have 
checked that the energy of the disordered low-temperature 
state of the system in our simulations is always slightly 
but significantly higher (by around 0.5% to 1%) than the 
energy of the fully ordered state if we include the full long- 
ranged interactions in calculating the energy 

The low-temperature states we see are governed mainly 
by the ice rule (though evidence of some local ordering 
of 4-clusters can be seen) and are probably macroscopic 
in number. This is why we observe an anomaly in the 
integrated entropy which we earlier attributed to a possi- 
ble ground state entropy. This low temperature "entropy" 
(as estimated from Figure 0) is around 10% lower than 
l/21og(3/2), which is itself an underestimate by around 
10%. Very much the same thing is likely to be true in 
the real system (dysprosium titanate) too: its true ground 
state is ordered but the system can almost never access this 
ground state. The measured ground state entropy here is 
closer to l/21og(3/2), in fact a bit more; it is probably 
a bit less than the true ground state entropy of nearest- 
neighbour spin ice, though. 

In the presence of a magnetic field, some interesting 
things happen. As reported earlier J3], some of the ob- 
served ground state entropy is recovered experimentally; 
we see this also in simulations (Figure 0). Since only 
the dipole moment couples to the field, the quantitative 
agreement in this curve is an additional confirmation of 
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Figure 8: Specific heat curves of dysprosium titanate (ex- 
perimental, and simulations give excellent agreement) and 
holmium titanate (experimental, and a simulation retain- 
ing dipole-dipole interactions to 12 nearest-neighbour dis- 
tances). By the time the ordering temperature is reached 
(which is expected to be the same for both systems) dys- 
prosium titanate is already frozen into an ice-like configu- 
ration from which it would find it hard to locate a lower- 
energy state. 



our model of co-existing dipole-dipole interactions and su- 
perexchange. The curves are similar in features and the 
amount of entropy recovered is also roughly the same. In 
stronger fields, sharp spike-like features start to show up in 
the experimental specific heat curves at low temperatures. 
Here, too, we find reasonable qualitative and quantitative 
agreement between simulations and experiment. All this 
confirms that our calculation of the dipole moment of the 
/ electrons and our supposition that the reduced energy 
scales are due to another interaction (superexchange) are 
correct, since the interaction with a magnetic field is purely 
magnetostatic. The experiments were done using powder 
samples, and the simulations show that the behaviour is 
strongly dependent on the direction of the field. To com- 
pare with powder averaged experimental results, we would 
need a very large number of simulations in random direc- 
tions, which we have not done to our satisfaction. 

The nature of the ground state also depends on the field 
direction, and with a sufficiently strong field and a suitable 
field direction the ground state may not even satisfy the ice 
rule. For instance, with a field along the z axis in Figure |3| 
(which corresponds to the [111] direction with the conven- 
tional unit cell) the ground state consists of three spins 
pointing into each upward tetrahedron and one pointing 
out of it. 

It is interesting to ask how the transition to a different 
ground state occurs as one slowly turns on a magnetic field 
at low temperatures. Figure nfi shows the result of doing 
this in a simulation at 0.2 Kelvin, for a field in the [ll\/2] 



"7 5 



^Rln2 
^•R(ln2-1/2ln3/2) 




T(K) 



Figure 9: The integrated specific heat per unit temper- 
ature in simulations without a magnetic field and with a 
half-tesla magnetic field. This shows the entropy gained 
over the ground state entropy. The entropy is expected 
to be i?log2 (dotted line) at high temperatures; the inte- 
grated value falls short of this, indicating a ground state 
entropy, but the ground state entropy is reduced in the 
presence of a magnetic field. Both the experimental data, 
taken from our earlier paper 0] , and the simulation results 
are plotted here for easy comparison. Based on simula- 
tions, we suggest that a magnetic field of around 3 Tesla 
should recover all or nearly all the ground state entropy. 



and [111] direction. The system seems to go through sev- 
eral magnetic transitions before reaching its fully polarised 
state. 

All these features would be averaged over in the experi- 
ments on the powder samples, and single crystals of these 
materials could turn out to be worth studying in their own 
right. 



4 An analogous six-vertex model 

If we examine the nature of states just above the transition 
in simulations of our holmium titanate model, we find that 
a large fraction of the tetrahedra are already in the ice- 
ruled state. This suggests that the sharp transition here is 
not the transition from paramagnetism to an ice-like phase, 
but a phase transition from a disordered ice-like state to 
an ordered or partially ordered state. 

To check that this is the case, we can try putting the ice 
constraint in by hand, and check that the phase transition 
is still reproduced at the same temperature. The model we 
get then is a form of six-vertex model on the "diamond" 
lattice, whose sites are the centres of the tetrahedra in the 
pyrochlore lattice. The six-vertex model has been widely 
studied on the square lattice; the diamond lattice, like the 
square lattice, has a coordination number of four and can 
be divided in two sublattices, but is three-dimensional. So 
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Figure 10: Specific heat in the presence of a 2 Tesla 
magnetic field, for various directions of the field. 



we study a system where one assigns arrows to bonds on 
the diamond lattice, such that each site (or "vertex") has 
two arrows pointing in at it and two pointing away: so six 
kinds of vertices are possible. But unlike conventional six- 
vertex models, we don't assign different weights to these 
six vertices, since all of them are really equivalent here; 
instead, the thermodynamics comes from interactions be- 
tween different vertices. In other words, we have a Hamil- 
tonian of the sort 



F = ^ = ^J(c(i),c(i),r i -r i ), 



(6) 



where c, is the configuration (a six- valued variable) of the i- 
th vertex, and J is the interaction energy of vertices i and j, 
which depends not only on their configurations but on the 
vector joining them (thanks to the underlying direction- 
dependent dipole-dipole interaction) . We have to calculate 
the pairwise J's appropriately. 

What we do is the following: we note that the sites on 
the diamond lattice fall into two sublattices, correspond- 
ing to up and down tetrahedra. First consider adjacent 
vertices (adjacent corner-sharing tetrahedra). The inter- 
nal interactions between these spins can be separated into 
nearest-neighbour interactions, which we can assume has 
already been taken into account via the ice rule, and next- 
neighbour interactions, which we can equally include by 
considering only next-neighbour vertices (that is, nearest- 
neighbour tetrahedra of like orientation) . We thus ignore 
interactions between nearest-neighbour vertices and con- 
sider interactions only between next-neighbour vertices, or 
nearest-neighbour vertices on a single sublattice. The in- 
teraction between two such vertices is the energy of inter- 
action between the two corresponding tetrahedra, as given 
by the sum of interaction energies of all pairs of spins. We 
ignore all further-range interactions. But the interactions 
already included, if carried out over all vertices over both 
sublattices, will actually double-count the pairwise spin- 



Figure 11: The growth of magnetisation in the simu- 
lation sample, for magnetic field in the [ll-\/2] and [111] 
directions, at a temperature of 0.2 K. 



spin interactions: we therefore also insert a factor of half. 
Formally, we can write 



H 



{id} 



E^ s 

{i,j} 



i , <Jj , i 



(7) 



and then break this sum up into nearest-neighbour terms, 
next-neighbour terms, and so on; throw out the nearest- 
neighbour terms because we have used them in enforcing 
the ice constraint; and group the next few terms into pairs 
where one spin belongs to one "up" tetrahedron, the other 
spin belongs to an adjacent "up" tetrahedron. 



H 



4 4 

/ J / J / J E(S am , S/3„, R Q — R/3 + x m — x„) 

{a,/3} 771=1 77 = 1 

-|-other terms (8) 



where the sum is over neighbouring "up" tetrahedra at 
sites a and /3, and we sum the interaction of each of the 
four spins at a with each of the four spins at 0. R Q is the 
position of spin 1 on tetrahedron a, likewise Rg, and x m 
are as in equation (0). All the "other terms" involve spins 
separated by three times the nearest-neighbour distance 
or more, so we drop them. We can do precisely the same 
grouping of terms for the "down" tetrahedra; we do both, 
and insert a factor of half. Thus our Ising spin Hamilto- 
nian is reduced to the vertex Hamiltonian ua) with appro- 
priately chosen interaction energies J, extending only to 
the nearest neighbour on the same sublattice. 

As before, we simulate this Hamiltonian. First, a word 
on how we do this. Flipping a single bond will not do: 
it will destroy the ice constraint on both adjoining ver- 
tices. We must find a closed loop, a set of bonds whose 
arrows lead from vertex to vertex and return to the start- 
ing vertex, and flip the whole loop at one go. Such "loop 
algorithms" have been discussed previously M and it has 



been pointed out that, in a six- vertex model, every line 
of arrows if followed must return to the starting vertex 
and every configuration is accessible via loop-flips alone, 
so a random loop-flip algorithm is ergodic. The earlier 
algorithms have several improvements and optimizations; 
however, they are concerned with conventional vertex mod- 
els where different vertices have different weights but don't 
interact, and cannot be completely translated to this sit- 
uation. We found it sufficient to merely pick up random 
starting sites, form loops randomly, calculate the energy 
difference, and flip them according to the Metropolis algo- 
rithm. 



The results are shown in Figure 12. It exhibits a phase 
transition at exactly the point where both the real sys- 
tem, and our model for it, do. This is a distinctly first 
order phase transition. Thus we have verified that the 
phenomenon driving this phase transition is not the for- 
mation of ice-like tetrahedra, but the further ordering of 
tetrahedra that have already attained ice-rule configura- 
tions; and we have displayed a fairly simple vertex model 
which shows the same features as our Ising pyrochlore. 

The ground state of this system would be expected to 
be fully ordered, but typically only partially ordered states 
are accessible. The argument is similar to that in the case 
of holmium titanatc, and the simulations bear this out. 
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lines in Figure G~3h. We ignore nearest-neighbour interac- 
tions for the same reason as earlier, ie that is taken care of 
by assigning an ice rule. 

The possible interactions are shown in Figure [IJ; for 
symmetry reasons, we need have only three interaction pa- 
rameters, all other non-zero interactions can be obtained 
from these by rotation, reflection, or inversion of one or 
both vertices (inverting a single vertex will simply change 
the sign of the interaction energy). If we calculate the 
interaction parameters from an assumed dipole-dipole in- 
teraction between magnets aligned along the edges con- 
necting the respective vertices, we obtain A = —8.6678, 
B = 10.753, C = -10.345 (arbitrary units). The ground 
state then looks as shown on the left of Figure 15; but 
a very slight alteration in the choice of A, B and C will 
give a ground state as shown on the right of Figure [15|, 
and the choice A — B = 2C may lead to rather interest- 
ing results. The model is probably worth studying both 
in its own right and because of the long historical inter- 
est square-lattice vertex models have held; but since it is 
not really related to the rest of this article, we postpone 
further discussion of it to a future work. 
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Figure 13: An ice model on the square lattice, with similar 
properties to the earlier, diamond-lattice model. Interac- 
tions are between diagonally-opposite vertices (grey lines). 



Figure 12: Comparison of specific heat curves for holmium 
titanate (the real system); our model of holmium titanate; 
the six-vertex model of section Ec and the six-state spin 
model of section j| 



5 A square lattice vertex model 

The sort of physics involved can perhaps be better seen in 
a square-lattice vertex model. Such models have been ex- 
tensively studied J8| , but the thermodynamics has typically 
arisen from assigning different weights to different vertices; 
instead we give the same weight to all vertices, but consider 
interactions between diagonally-opposite vertices (the grey 
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Figure 14: The three possible interactions between neigh- 
bouring vertices. All other possibilities are symmetry re- 
lated, or zero. In particular, reversing all arrows on a single 
vertex will simply change the sign of the interaction. 
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Figure 15: With interaction energies A, B, C between 
vertices calculated from dipole-dipole interactions between 
spins aligned along their edges, we obtain a ground state 
as on the left. But a slightly different choice of weights 
will yield a ground state as on the right, and there is the 
possibility of a "level crossing" between the ground states 
for a choice A — B = 2C ' . 

6 Multistate spin model on an fee 
lattice 

In our vertex model earlier, we had two sublattices, and in- 
teractions only within a single sublattice. Apart from the 
ice-rule constraint, the two sublattice could just as well 
be noninteracting. So the next logical step is to separate 
the two sublattices. We consider an fee lattice, with a six- 
valued variable at each site. Only nearest-neighbour inter- 
actions are considered, and as before, the value of the inter- 
action is determined from the underlying pyrochlorc Ising 
variables. The major difference with the vertex model case 
is that we have now forgotten about the "down" tetrahe- 
dra: an arbitrary configuration of four neighbouring "up" 
tetrahedra would not satisfy the ice rule for the enclosed 
"down" tetrahedron, but we are no longer worrying about 
that now. 

It turns out that the dynamics of interaction between 
these "up" tetrahedra takes care of that for us. The sys- 
tem displays a phase transition at very nearly the same 
temperature as the vertex model and the Ising pyrochlore 
(Figure 12), and at temperatures just above the transition 
the configuration is such that nearly all the "down" tetra- 
hedra in the correspondingly-configured pyrochlore would 
satisfy the ice rule; at zero temperature the system is par- 
tially ordered, in exactly the way we observed in holmium 
titanate, but the partially ordered states in this case really 
are the ground states. 

The partial ordering is now more easily visualised with 
the conventional cubic unit cell rather than the paral- 
lelepiped which we used. Note first that each ice-ruled 
state of a tetrahedron has a dipole moment, perpendicular 
to the side connecting the two inward-pointing spins and to 
the side connecting the two outward-pointing spins. If one 
looks at the cubic unit cell, we can see that the six allowed 
values of this dipole moment are along the three edges of 
this cube: so what we have is a six-state magnetic model 
on an fee lattice where each spin can point along one of 
the Cartesian axes. In the ground state, one of these axes 
is picked out, so that each spin points along the same line 



in one of two opposite directions; each plane perpendicu- 
lar to this direction is antiferromagnetically ordered; and 
perpendicular to this plane, the ordering is random. 

It is tempting to use the total dipole moments of the 
tetrahedra as the site variables, and for the interaction 
simply to use their mutual magnetic interactions, since we 
know that the dipole-dipole interaction favours antiferro- 
magnetic ordering in planes perpendicular to the spins; but 
it turns out that this is not the ground state of such a sys- 
tem. Only by using the actual interaction energies of the 
tetrahedra do we obtain such a ground state. However, 
there is an obvious connection between this system and a 
well studied-problem, which we turn to in the next section. 

7 Ising model on the fee lattice 

The nearest-neighbour antiferromagnetic Ising model on 
the fee geometry has been studied by several authors, and 
its ground state is known to have exactly the sort of or- 
dering we are considering. 

The ordering is easy to understand if there is a bit of 
anisotropy in the system: consider Figure 16L where we 
have an fee crystal, and within the x-y plane and planes 
parallel to it there is an antiferromagnetic interaction J be- 
tween nearest neighbours, but out of the plane there is an 
interaction J' < J. Then the planes prefer to order antifer- 
romagnetically, but adjacent planes have zero interaction 
energy regardless of their relative ordering, so the ordering 
along the z axis is random. The interesting thing is that 
this remains true even when J' = J: the only change is 
a new factor of 3 in the degeneracy because of the new 
rotational symmetry of the system. 




Figure 16: Ising antiferromagnet on an fee lattice, with 
nearest neighbour interaction J (solid grey lines) in the 
plane and J' < J (dotted lines) between planes. Each 
plane orders antiferromagnetically but they stack up in a 
random manner. Surprisingly, this is also the ground state 
when J' = J. 

The order of the transition is also of interest. In the 
isotropic case, it is a first order transition. With strong 
anisotropy (a weak inter-plane coupling), however, we 
would expect a second-order transition because the sys- 
tem effectively is like weakly interacting 2D Ising systems, 



which have a second order transition at the Onsager tem- 
perature. In fact, simulations suggest that for J' close to 
J, the transition is first order, but for J' somewhat less it 
is second order, and for J' = 0.6J the transition temper- 
ature is almost exactly the Onsager temperature. Since 
the planes have a zero interaction energy in the ground 
state and a very small interaction energy at low tempra- 
tures, they behave like almost uncoupled 2D Ising systems. 
Thus this system seems to exhibit either a first order tran- 
sition or a second order transition with the same sort of 
ground state, depending on the parameters. 



8 Conclusion 

We have clarified the true nature of the ground states of 
the Ising pyrochlores holmium titanate and dysprosium ti- 
tanate. We have pointed out that the ice-like behaviour 
of dysprosium titanate seems to arise not from a macro- 
scopic degeneracy of the true ground state, but from its 
inaccessibility in practice, and consequently the tendency 
of the system to fall into one of a large number of slightly 
excited ice-like states. In holmium titanate, the order- 
ing temperature is higher than the expected ice-formation 
temperature; here, too, the system gets stuck into excited 
states, but these are partially ordered states and the model 
system shows a clear phase transition. By rigidly enforc- 
ing the ice constraint, we show that this transition exists 
independently of the broad ice-state crossover in spin ice, 
and we exhibit several models, including the well-known 
fee Ising model and a diamond lattice vertex model, which 
undergo a similar phase transition. Analogous to this ver- 
tex model we also exhibit a square lattice vertex model 
which has differently ordered ground states depending on 
what interaction parameters we choose, and which we hope 
to examine further sometime in the future. 

In addition, we have looked at what happens to dyspro- 
sium titanate when a magnetic field is applied, and com- 
pared our conclusions to available experimental data; we 
have reproduced earlier experimental data for a weak field 
and see similarities in our simulations and the strong-field 
data, further confirming our underlying model. We see 
strongly anisotropic behaviour in these systems and pre- 
dict some interesting results for possible experiments on 
single-crystal samples. 

This work appears as part of the Ph. D. thesis of R. S.M 
Note. After this manuscript was prepared, we learned 
that the ordering discussed by us has recently been ob- 
served in simulations by den Hertog et aZ. [l0| 
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